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In a model of self-organized criticality unstable sites discharge to just one of their neighbors. For 
constant discharge ratio a and for a certain range of values of the input energy, avalanches are simple 
branchless Polya random walks, and their scaling properties can be derived exactly. If a fluctuates 
widely enough, avalanches become branched, due to multiple discharges, and behave like those of 
the stochastic sandpile. At the threshold for branched behaviour, peculiar scaling and anomalous 
diffusive transport are observed. 
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An avalanche is a cascade of a large number of micro- 
scopic events. Generally it is triggered by a single event 
at a point. This creates similar events in the neighbour- 
hood which activate the neighbours at further distances. 
Thus, an avalanche is a simple example of branched 
growth process, of which several physical examples exist 
0. Particularly interesting are situations in which the 
avalanche sizes have a scale free distribution, typically a 
simple power law, signifying the existence of long-range 
correlations as in critical phenomena. This is the case of 
avalanches in the phenomenon of Self-Organized Critical- 
ity (SOC) H where long-ranged spatio-temporal corre- 
lations spontaneously emerge in non-equilibrium steady 
states of slowly driven systems with nonlinear local re- 
laxation mechanisms J^Q]. 

Sandpiles are prototype models of SOC. rij grains re- 
side at the i-th site of a regular lattice. The role of the 
external drive is to trigger transport processes through 
the system by adding single grains of sand at a time at 
randomly selected sites: ni — > rij + 1. If > n c the sand 
column at i topples and grains are distributed to the 
neighbourhood. In the BTW model, the grain distribu- 
tion process is deterministic since each neighbouring site 
gets one grain In the stochastic two-state sandpile 
model grains are transferred to randomly chosen neigh- 
bouring sites In spite of their similarities, the BTW 
model and the stochastic sandpile are now believed to 
belong to different universality classes |J. Indeed, very 
recent studies show that the BTW model has a multifrac- 
tal behaviour whereas standard finite size scaling 

works well for the two-state stochastic sandpile [p|~|TT| . 

Characterizing features of all these sandpile models are 
(i) the normal diffusive dynamics of the particles and (ii) 
the branching of the toppling process. A grain moves a 
unit distance in a toppling and its resultant motion un- 
der different topplings in different avalanches is diffusive. 
This implies that the average number of topplings in an 
avalanche grows as a quadratic power of the system size, 
L, in all isotropic sandpile models. Secondly in all mod- 
els, the toppling condition is made in such a way that a 



single toppling can excite more than one neighbors, which 
ensures that an avalanche is a branched process. Normal 
diffusive transport and branching are believed to be very 
basic ingredients of SOC. 

In this paper, by studying an original energy activa- 
tion model with stochastic discharge mechanism, we find 
that neither of these ingredients (i) and (ii) as stated 
above, is necessary for SOC behaviour. Indeed, we show 
that avalanches in SOC can be linear as branchless ran- 
dom walks. In this case their scaling behaviour is dif- 
ferent from that of branched avalanches, Eulerian ran- 
dom walk models are also linear avalanche models stud- 
ied before . and the relative probability distributions 
do not become independent of L, as this tends to in- 
finity. By allowing the discharge ratio to fluctuate in 
a progressively wider interval, we are also able to trig- 
ger a branched avalanche behaviour, which falls in the 
universality class of the stochastic sandpile || . Right at 
threshold for this branched behaviour the diffusive trans- 
port becomes anomalous, i.e. the average avalanche size 
grows as a power higher than 2 of L. 

Our model has two parameters: the input amount S 
and the discharge ratio a. An amount of energy ej, which 
can vary continuously, is accumulated at each site i of a 
square lattice box of size L. The system is driven by 
injecting an amount of energy S to a randomly selected 
site. Every lattice site has a limiting capacity e c (= 1) for 
the maximum amount of energy storage. If at any site 
the energy > e c , the site i activates and undergoes a 
relaxation process. In a relaxation a fraction a of the site 
energy is discharged to only one of the neighbouring sites, 
j, which is selected randomly: tj — > tj + ae^. The rest of 
the energy remains in the relaxing site: ti — * (1 — a)ej. 
If the energy at the receiving site j exceeds the thresh- 
old, that site also relaxes. The possibility of multiple 
relaxations arises when site i remains active even after 
relaxing. In such a case further relaxations occur in the 
subsequent stages. Upon relaxing, an active site at the 
boundary may drop a fraction a of its energy outside the 
system. This ensures that the system reaches unique 
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FIG. 1. (a) A steady state energy configuration in the LAM (a = 1/3, 8 = 1/5) within a circular region on a square lattice 
of size L = 81. (b) Plot of the correlation function C(r) for the same LAM as in (a) but for L — 257. 



stationarity after repeated driving. Our simulations 
show that this is the case, independent of starting con- 
figuration as required in SOC, for any particular selection 
of (a, 8). We note that this model is a special case of the 
class of models called the abelian distributed processors 
model @. 

First we consider a constant value of a. In the steady 
state, 1 — a < e, < 1 since each site relaxes at least once. 
This implies that a site must activate after receiving at 
least energy a due to an activation at the neighbouring 
site. Energy activation events then take place one after 
another like in a step sequence of a Polya random walk 
(RW) ju). Once started, such a RW must terminate at 
the boundary, after dropping some energy outside the 
system. We call this the linear avalanche model (LAM). 

When 8 < a there must be some sites which do not ac- 
tivate upon receiving the input energy, since a RW drops 
at least a energy outside the system. To eliminate corner 
effects we use a circular region, within a square box of 
size L (odd), of radius R = (L — l)/2. Fig. 1(a) shows 
an energy configuration in the steady state where the 
sites having energy larger than the average, (e), are rep- 
resented by black dots. Correlated regions of such sites 
are observed as connected clusters of black dots. The cor- 
relation function C(r) is defined as the probability that 
a black dot at a distance r from the centre belongs to the 
same cluster containing a black dot at the centre. We ob- 
serve a power law decay C(r) ~ with £ = 0.37 ± 0.02 
which is the signature of the critical correlation devel- 
oped in the steady state (Fig. 1(b)). However a site 
energy correlation function like: (e(O)e(r)) — (e 2 ) aver- 
aged over all sites has an almost uniform small negative 
value except when r ~ R. 

The probability density of the site energies D(e,R) in 
the steady state is very well fitted with a generalized 
Lorentzian peak in the range 1 — a < e < 1 around its 
average (e(. D follows a scaling with R of the form: 



D{e)lR = a*/{{e-(e)YR* + b 2 Y (1) 

where, for example, a rs 0.195,6 w 0.226 and c ks 1.33, 
with a = 1/3 and 8 = 1/5 (Fig. 2). From Fig. 2(a) we 
see that the site energy distribution is symmetric about 
its average value at (e) and the most probable value co- 
incides with the average value. Assuming that the e > 
is less than 1 — 8, the fraction of sites which donot ac- 
tivate on receiving the input energies are greater than 
those which activate. This situation cannot be stable 
since these low energy (e < 1 — 8) sites absorb energy 
at a greater rate from the external drive than the high 
energy sites (e > 1 — 8) and the average energy will push 
up. Similarly if (e) is greater than 1 — 8 there would be 
faster dissipation of energy through the boundary which 
will push down the average energy. Therefore it is likely 
that the stable states correspond to (e) = 1 — 8 where 
the rates of absorption by the low energy sites and the 
rate of dissipation by the inputs at the high energy sites 
are equal. Our numerical results strongly supports this 
result: 

{e(a,8)) = l-8 for 8 < a (2) 

For a < 8 < 8 c (a) every RW drops exactly an amount 
of energy 8 outside the boundary. D(e, R) is a delta func- 
tion at its average value (e) such that a site after activa- 
tion retains the same energy which it had before receiving 
the external input 8. This implies 

(l-a)((e)+8) = (e) (3) 

which gives (e(a,8)) = (1/a — 1)8. The average en- 
ergy increases to 1 when the input energy is increased 
to 8 c (a) — a/(l — a). This is the limiting situation for 
the branchless avalanches. Beyond this limit multiple 
discharges start and avalanches cease to be branchless. 
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FIG. 2. D(e) for (a) LAM and (b) BAM. The continuous curves are fits to a generalized Lorentzian function in (a) and to 
a power law in (b). 



The size of the avalanche is generally measured by 
the total number of relaxations (s), the life-time of the 
avalanche (t) and the radius of the avalanche (r). Using 
{x : s, i, r} we assume x ~ x' lx '<* , where the j rx is the cut- 
off exponent for the measure x, and the usual finite size 
scaling forms for the probability distribution functions of 
avalanche sizes: Prob(x, L) ~ L~" x f x ( L * rx ) , where, 

fx(y) ~ y~ Tx in tnc umit 01 2/ -* ° s ives T ^ = Px/irx, 

in case Prob does not maintain an L dependence for 
s « L. The 7 exponents are connected by the rela- 
tions: j x ' x = 7x'x"7x"x and ~/ xx > = (t x > - 1)/(t x - 1). 

The steady state of the LAM is related to the first 
passage problem of RW's. The avalanche size and the 
life-time are the same as the number of steps s taken by 
the walker before dropping through the boundary. The 
distribution of the avalanche sizes can be calculated in 
the following way. The probability to start a RW within 
a distance AR from the boundary is 2itRAR/ttR 2 ~ 
AR/R. A RW moves a distance ~ s 1 / 2 in s steps. There- 
fore, the probability that an arbitrary walker makes s 
steps or less before reaching the boundary is the prob- 
ability that an arbitrary site is within a distance s 1 / 2 
from the boundary, which is ~ s 1 / 2 / R. The probabil- 
ity that an arbitrary walker makes precisely s steps to 
reach the boundary is s -1 / 2 / R. Thus, in terms of L, 
Prob(s,L) ~ s- Ts /L for s « L with r s = 1/2. Our 
simulations show that upon increasing the system size 
the effective exponent t s (L) gradually decreases to its 
asymptotic limit of 0.50(1). Since the avalanches are ran- 
dom walks, the cut-off for s must be proportional to L 2 , 
and therefore 7 rs = 2. Normalization of Prob(s, L) gives 
(3 S = 2. The relation r s = {3 s /"/ rs is not satisfied by the 
LAM avalanches, because Prob maintains an L depen- 
dence for s << L. This unusual dependence is due to 
the fact that LAM avalanches must necessarily reach the 
boundary in order to extinguish. As a rule, this is not 
the case for other SOC systems. 



We compare these results with two cases studied in the 
literature. A model of SOC where avalanches are branch- 
less walks is the Eulerian walkers model. In this model 
each site of a regular lattice has an outgoing direction, 
which is one member of a set of outgoing bonds asso- 
ciated with this site. The walker leaves the site along 
the outgoing direction but changes the outgoing direc- 
tion sequentially to the next member in the set of out- 
going bonds 0|. In the Euler walker case, the depen- 
dence of Prob(s, L) is of the form L~ 2 f(sL~ 2 ) just sim- 
ilar to our LAM where s is the number of steps to the 
boundary. Secondly for the random walks on a square 
with absorbing boundary, the exact scaling function for 
the distribution of number of steps to the boundary is 
P(s) ~ L~ 2 f(sL~ 2 ) has been calculated. It has been 
shown that the scaling function f(x) can be expressed 
explicitly in terms of the Jacobi theta function . 

Next we study a situation when the discharge ratio a 
is a random variable and a fresh value for it is drawn 
from a uniform random distribution in {0, 1} each time 
a site relaxes. Since now the energy released after acti- 
vation may be arbitrarily small, multiple relaxations are 
quite frequent, and this leads to branching. We call this 
the branched avalanche model (BAM). The probability 
density D(e) of the site energies in the steady state has 
now little dependence on S as well as on L. D(e) is fitted 
best to a form D(e) = D a + D^, with D a = 0.41(1), 
D 1 = 1.59(1) and n = 1.70(2) (Fig. 2). The average 
energy per site (e) is however observed to have an L de- 
pendence as: (e) = f M — CL~ 3 / 4 with — 0.6365(3) 
and C = 0.21(2). 

The probability distributions for avalanche sizes and 
life-times follow the finite size scaling form of Prob very 
well. The data for Prob(s, L) for three different system 
sizes L = 129, 513 and 2049 collapse well for 7 r;s = 2.75 
and P s = 3.6. This gives a value for r s = 1.31. Similarly, 
we obtain r t = 1.51(3) and j rt = 1.50(3). The average 
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size (s(L)) ~ L Va with v s = 2 and similarly = 0.76(3) 
are obtained. The values of these exponents are very 
close to those of the two-state model, indicating that the 
BAM may well belong to the universality class of the two- 
state model ||. In the BAM, avalanches can extinguish 
within the boundary. 



(e)(l + W2 + w% + wl + + W2) + W2<5 which in the s — > oo 

limit gives (e)/(l— K/2). If this site now releases a fraction 
Wi of its energy, the amount left is (e) ((1 — Wi)/(1 — W2)) 
which has to be greater than one for a second relaxation 
to take place. Now since (e) can be at most (1 —6), one 
must have w c > 5/(2 — 5). 
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FIG. 3. Variation of the average branching number 
r(w,L) with w. The curves become steeper as the system 
size L increases from 33 to 512. The threshold for BAM is at 
w c = 0.125 ±0.010. 

The LAM and the BAM can be regarded as the two 
extremes of a range of situations characterized by a pro- 
gressively increasing width of the interval in which a is 
allowed to fluctuate. We introduce a parameter w as the 
size of the window for the stochastic a, which is now cho- 
sen randomly within the range {1/2 — w/2, 1/2 + w/2} 
with uniform probability. When w = we have the LAM 
with a constant a = 1/2, whereas when w = 1 we have 
the BAM. A full characterization of all dynamical regimes 
for < w < 1 turns out to be extremely challenging. 
However, starting from w = 1 and progressively lowering 
this parameter, we are able to identify a whole region 
w c < w < 1 in which BAM behaviour holds, and to 
characterize peculiar, novel scaling at the threshold for 
BAM behaviour, w — w c . In the interval < w < w c for 
sure LAM behaviour prevails in a whole neighbourhood 
of w = 0. However, a sort of double degeneracy of the 
stationary state arises for higher w's, making the scaling 
analysis quite difficult. This degeneracy is an interesting 
phenomenon in itself, worth further investigations, and 
indicates that the LAM-BAM transition is a very com- 
plex process. 

A lower bound for w c can be estimated by calculat- 
ing the maximum amount of energy a site can receive 
which is needed for at least two relaxations. Suppose 
all sites have the same energy (e) and W\ = 1/2 — w/2 
and W2 — 1/2 + w/2. Then on adding an amount 6 
of energy at a site, the maximum amount of energy 
with which a site relaxes after the s-th step of a ran- 
dom walk is: ....(w 2 (w2(w2(6 + (e)) + (e)) + (e)))..)... = 
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FIG. 4. The data collapse of the difference in the average 
size of the avalanche (s(w,L)} and that in the BAM (s(l,L)) 
is plotted for L = 65, 129 and 257. The inset shows the vari- 
ation of {s(w,L)) ~ L"'^ where v s {w) = 2.75, 2.01, 2.00 for 
w = w c , and 1 respectively. 

The basic difference in the limiting cases w = and 
w = 1 is the branching of the avalanches. Therefore, a 
quantity which one can monitor in order to describe the 
passage from LAM to BAM behaviour is the ratio (s)/(t). 
If on average an avalanche has size (s) and life-time (t), 
(s) / (t) measures the number of branches, once assumed 
that each branch has a duration (t). Therefore, we define 
the scaled branch number as the order parameter: 



r(w,L) = ((s)/(t)-l)/L^-^ 



(4) 



and define it as an order parameter in the limit L — *• 00. 
By definition this order parameter must be finite and 
nonzero for w > w c . It also turns out to be equal to 
zero for w < w c . Fig. 3 shows that for any L, r(w,L) 
is very close to zero for w < w c , but it increases very 
fast for w > w c . Upon making L larger, r(w,L) be- 
comes progressively smaller for w < w c , whereas it rises 
at a faster rate while w > w c . Using 5 = 1/8 we iden- 
tify w c — 0.125 ± 0.005. The sharp increase of r(w,L) 
as w — > u>+ is well fitted as r(w,L) ~ (w — w c )@ , with 
[3 = 1.2 ±0.1. 

The average size of the avalanche (s(w,L)) has a 
sharp peak at w c (L) and it depends on L as: w c (L) = 
w c + 0.39L~ 1/2 with w c = 0.123. The data collapse 
well when ((s(w,L)) — (s(l,L)))L~ 5 is plotted versus 
(w - w c )L - 34 . This indicates that (s(w c ,L)) - L 2 - 75 , 
which implies an anomalous diffusive transport at thresh- 
old (Fig. 4). Such anomalous transport has never been 
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reported before for a SOC model, to our knowledge. At 
w = w c the branch number (s) / (t) grows as L . So, the 
threshold regime with anomalous diffusive transport is 
still characterized by infinite branching of the avalanches. 

The dependence of {s(w c , L)) ~ L 2 - 75 at w — w c seems 
quite surprizing. We provide the following tentative ex- 
planation for the same: At any nonzero value of w some 
avalanches are purely random walks (i.e. s — t) and oth- 
ers are branched. If fnw and Jba = 1 — Jew are the frac- 
tions of avalanches which are linear and branched, then 
for w = 0, fnw = l and it decreases as w increases, conse- 
quently Jba increases from its zero value. Both fractions 
become very close to 1/2 at around w = 0.4 and beyond 
this value they vary slowly to fuw ~ 0.6 and fsA ~ 0.4 
at w = 1. We also observed the fraction of RWs which 
terminate within the boundary of the system and donot 
drop out any energy outside the system. We find that the 
fraction of such random walk avalanches increases from 
zero very sharply to almost 0.9 at around u>=0.15. This 
implies that around w — w c most of the 8 input energies 
corresponding to those RWs which terminate within the 
boundary get stored in the system. Therefore, it is the 
branched avalanches which take out this extra stored en- 
ergy from the system - consequently their sizes are bigger 
and perhaps for this reason the average size varies as a 
larger power of L like 2.75. 

At w — w c , the avalanche size distribution D(s) vs. 
s has two regions with two characteristic sizes and 
s 2 . While s 2 ~ L 2,75 is the usual cut-off size for the 
avalanches, s\ ~ L 2 is an intermediate size. For s < si, 
tI = 0.64 whereas for the second region s). < s < s 2 the 
value of t 2 is around 1.45. For w < w c but near to it, 
s\ r~j L 2 whereas s 2 grows as w approaches w c . How- 
ever, for w > w c , s\ decreases as s\ ~ (w — w c )~ 2 ' 25 as 
the deviation (w — w c ) increases and r 2 also continuously 
decreases to 1.3, while tI — > as L — > oo. 

To summarize, in a stochastic energy activation model 
an active site transfers the energy to only one randomly 
chosen nearest neighbour. The avalanches are linear 
when the transfer amount is narrow distributed and are 
branched when the transfer is broadly distributed. While 
the exponents of the toppling distribution of the branch- 
less avalanches can be exactly determined, the universal- 
ity class of the branched ones appears compatible with 
that of the stochastic two-state sandpile. A transition 



between the two regimes is observed by tuning the size 
of the window of the stochastic transfer ratio. At the 
transition point, the diffusion mechanism induced by 
avalanches is anomalous. 
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